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ABSTRACT 

It has been known for over a decade that many four-image gravitational lenses exhibit 
anomalous radio flux ratios. These anomalies can be explained by adding a clumpy 
cold dark matter (CDM) component to the background galactic potential of the lens. 
As an alternative, Evans & Witt (2003) recently suggested that smooth multipole 
perturbations provide a reasonable alternative to CDM substructure in some but not 
all cases. We generalize their method in two ways so as to determine whether multipole 
models can explain highly anomalous systems. We carry the multipole expansion to 
higher order, and also include external tidal shear as a free parameter. Fitting for 
the shear proves crucial to finding a physical (positive-definite density) model. For 
B1422-I-231, working to order /cmax = 5 (and including shear) yields a model that is 
physical but implausible. Going to higher order (fcmax ^ 9) reduces global departures 
from ellipticity, but at the cost of introducing small scale wiggles in proximity to 
the bright images. These localized undulations are more pronounced in B2045-f265, 
where fcmax ~ 17 multipoles are required to smooth out large scale deviations from 
elliptical symmetry. Such modes surely cannot be taken at face value; they must 
indicate that the models are trying to reproduce some other sort of structure. Our 
formalism naturally finds models that fit the data exactly, but we use B0712-I-472 
to show that measurement uncertainties have little effect on our results. Finally, we 
consider the system B1933-I-503, where two sources are lensed by the same foreground 
galaxy. The additional constraints provided by the images of the second source render 
the multipole model unphysical. We conclude that external shear must be taken into 
account to obtain plausible models, and that a purely smooth angular structure for 
the lens galaxy does not provide a viable alternative to the prevailing CDM clump 
hypothesis. 

Key words: galaxies: haloes — galaxies: structure — dark matter — gravitational 
lensing 



1 INTRODUCTION 

It has been suspected since tlie time of Newton that a light 
ray would be deflected when it passes near a massive ob- 
ject. This phenomenon, now known as gravitational lensing, 
provided one of the earliest and best tests of the validity of 
general relativity. Since that time, lens theory has matured 
into an active field of astronomy (see Kochanek et al. 2004, 
for a recent review) . Probing the structure of galaxies is one 
of the many areas to which lensing has been applied over the 
past fifteen years (e.g., Kochanek 1991; Keeton et al. 1998). 
Because lensing is sensitive to all mass, it is possible to study 
both dark and luminous components within galaxies. 

In the case of strong lensing where the light source is 
multiply imaged, we can use the positions and fluxes of 
the images to study small-scale structure within the lens 



galaxy. This technique can be most readily applied to four- 
image systems in a cusp or fold conflguration, which occurs 
when the angular separation between the lens and source 
is small. For a cusp lens we expect the three brightest im- 
ages to satisfy the magniflcation sum rule (e.g., Schneider & 
Weiss 1992; Mao & Schneider 1998; Keeton et al. 2003), 

Ml + + ~ 0. (1) 

In the fold case an analogous relation holds for the two 
brightest images (Gaudi & Fetters 2002; Keeton et al. 2005): 

^ll+^l2~0■ (2) 

Although these relations should hold for all smooth lens 
potentials, there are a number of observed systems for which 
they are violated (Mao & Schneider 1998; Keeton et al. 2003, 
2005). Since the magnifications are not directly observable, 
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we refer to systems tliat violate tlie sum rule as exhibiting 
flux ratio anomalies, with the observed fluxes being related 
to the magnifications by the flux of the source. 

One might argue that violations of the magnification 
sum rule originate in electromagnetic effects of the interstel- 
lar medium on light emitted by the lensed source. However, 
dust extinction is negligible at radio wavelengths, which are 
much larger than the typical size of a dust grain. The lack of 
wavelength dependence in radio flux ratios rules out electro- 
magnetic effects as an explanation for radio anomalies (see 
e.g., Kochanek & Dalai 2004; Keeton et al. 2003, 2005, and 
references therein). 

Explaining the observed anomalies therefore requires 
the addition of small-scale structure to the gravitational po- 
tential of the galaxy (Mao & Schneider 1998). A possible 
candidate for this substructure emerged through the work 
of Moore et al. (1999) and Klypin et al. (1999), whose nu- 
merical simulations predicted a quantity of cold dark matter 
(CDM) halos with masses ~ lO'' — 10** Mq , greatly exceeding 
the observed numbers of such objects. This result motivated 
Metcalf & Madau (2001) to consider how the abundance 
of predicted CDM substructure might affect lensing. They 
pointed out that even if the small halos arc completely dark 
— invisible to standard luminosity-based studies — they 
could still affect lens flux ratios and perhaps explain the ob- 
served anomalies. Indeed, DalaJ & Kochanek (2002) found 
that putting ~ 2% of the mass in ~ W^Mq halos could re- 
produce the observed flux ratios for seven anomalous lens 
systems, while broadly matching the predictions of the nu- 
merical simulations. But does this mean that the flux ratio 
problem is really solved, and that the "missing" CDM satel- 
lites have been found? Are there other plausible models that 
can solve the flux ratio problem? 

Possibilities may include stellar microlensing, and com- 
plex structure such as isophote twists or triaxiality in the 
lens galaxy. Since the radio-emitting regions of the QSO 
sources we will consider have a much larger angular scale 
than the Einstein radius of a typical star in the lens galaxy, 
microlensing can be eliminated as a potential cause of flux 
ratio anomalies. That leaves the question of whether models 
that alter the global structure of the lens potential offer a 
viable explanation of flux ratio anomalies. Our goal is to see 
whether we can fit four anomalous radio lenses using models 
with general but reasonable angular structure. 

We begin with the self-similar multipole model of Evans 
& Witt (2003). In this framework, the potential of the lens 
galaxy is described by a generalized isothermal model whose 
angular dependence is expressed as a Fourier series. The 
multipole coefficients are determined by fitting the observed 
image positions and flux ratios. The truncation order of the 
series is chosen so that the matrix of constraints is square. 
Evans & Witt (2003) found that such a model could ex- 
plain some but not all observed lenses. In particular, for 
B1422+231 (the only radio anomaly they studied), they 
could find a physically acceptable model only by inflating 
the errorbars on the data, and even that model was highly 
non-elliptical and implausible. 

In this paper we extend the multipole formalism to in- 
clude external shear tidal distortions from objects near 
the lens galaxy (e.g., Keeton et al. 1997) — and higher or- 
der multipole modes. Shear in particular will prove essential 
for obtaining sensible solutions. While this is not the most 



general method, it is mathematically simple and may be a 
reasonable alternative to substructure in some lenses. 

As a test case, we first apply our model to Q2237+0305, 
which is not anomalous at radio wavelengths. Then, we ap- 
ply the model to B1422+231, B2045+265, B0712+472, and 
B1933-(-503, which are all highly anomalous. An exhaus- 
tive study of the known radio anomalies would also include 
B 1555-1-375, but the position of the lens galaxy (a key in- 
gredient in our formalism) is unknown. 



2 METHODS 

We begin by writing the convergence (dimensionless surface 
mass density), k, and lens potential, -ih, for a galaxy with 
a flat rotation curve and arbitrary angular structure. This 
model is often referred to as a generalized isothermal model 
(Witt et al. 2000; Zhao & Pronk 2001; Evans & Witt 2001, 
2003): 

K{r,e) = ^G{e); i^{r,e)=rF{0). (3) 
Noting that V^V — 2k, we find that F and G are related by 

G{e) = Fie) + F"ie). (4) 

For a given source position, u = {u, v), we can find the image 
positions, x = {r cos9,r sinO), via the lens equation, 

u = X - VV'(x) . (5) 

The inverse magnification of an image at x is given by 

An important property of the lens potential is the criti- 
cal curve, along which the magnification is formally infinite. 
The critical curve in the image plane maps to the caustic 
in the source plane, which marks the transition between 2- 
image and 4-image systems. We see from equation (6) that 
the critical curve is given by rcrit(^) = G{6), which is equiv- 
alent to the isodensity contour with k = 1/2. 

2.1 Multipole Lens Model of Evans &c Witt (2003) 

Let us write the angular part of the potential, F, as a mul- 
tipole (Fourier) expansion in 6, i.e. 

femax 

F{e) s y -f ^ (ofc cos kO + bk sin kO) , (7) 

fc=i 

for some appropriate fcmax. The Fourier coefficient ao is re- 
lated to the Einstein radius by _REin = ao/2. We can find the 
unknown source position, u, and the Fourier coefficients, 
and bk, by introducing observational constraints, viz. the 
image positions and fiux ratios. 

From equations (3) and (7), the lens equation (5) be- 
comes 

fcmax 

u = n cosOi - Y<^oiei) - 2 [akakiOi) + bkPkiOi)] (8) 

k=l 

fcmax 

V = ri smOi - yao(6li) - £ [afc«fc(^0 + bk$k{Oi)] (9) 

k=l 
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where Z = 1, 2, . . . , n is the image number. The functions 
Ok, Ok, Pk, and $k are defined in Evans & Witt (2003), 
equations (13) and (14): 



Qfc(^) ~ cos 6 cos k9 + k sin 6 sin k6 , 

Qfe {0} = sin 9 cos kO — k cos 9 sin k9 , 

/3k (9) = cos 9 sin kO — k sin 9 cos k9 , 

$k {0) = sin 9 sin k9 + k cos cos kO . 



(10) 

(11) 
(12) 
(13) 



Another set of constraints comes from the flux ratios. 
Relative to image n, we have 



fni 



tin 
1^1 ' 



(14) 



We then use (4), (6) and (7) to obtain 



ifni - l)r„r, = y7o(ft) + Yl [«fc7fc(^0 + bk5k{ei)] , (15) 

fc=i 

where / = 1,2, ...,n— 1. The functions 7fe and Sk are defined 
by Evans & Witt (2003), equation (18): 

IkiOi) = {1 -k'^)[fniri cos k9n-r„ cos k9i] (16) 
5k{0i) = {1 -k^)[fniri sin k9n -rn sin k9i]. (17) 

Wc can combine equations (8), (9) and (15) into a single 

matrix equation: 



(18) 



where x is the (2femax + l)-dimensional vector 
of parameters for which we are solving; x = 
(u, w, ao, a2, 62, • • • , ^fcmox)- Wc drop ai and 61, 

which represent an unobservable translation of coordinates 
in the source plane. The (3n — l)-dimensional vector b 
contains the observed image positions and flux ratios; 



= {xi,...,x„,yi,...,y„, 

ifni - l)r„ri, . . . , (/n,n-i - l)r„r„_i), 



(19) 



where xi = n cos 9i and yi = ri sin 9i . The (3n — 1) x (2fcmax + 
1) matrix. A, is defined in equation (22) of Evans & Witt 
(2003): 



A = 



1 

1 
1 



1 




aoi 

aon 
&01 

701 



021 
0121 

721 



1321 

/32n 
P2I 



I32n 

hi 



(20) 



. 7o,n-l 72,n-l <52,n-l 



where aki = ak(9i) etc. Evans & Witt (2003) choose fcmax 
such that A is squaxe. We then have fcmax = 5 for a 4-image 
system (n = 4) . With this choice of fcmax we can simply mul- 
tiply equation (18) by A~^ to solve for x, provided that A is 
non-singular. To ensure numerical stability, however, Evans 
& Witt (2003) solve for x using singular value decomposition 
(SVD). 



2.2 The Minimum Wiggle Model 

There are two main limitations of the method of Evans & 
Witt (2003). On a technical point, their requirement that 
A be square prevents one from probing the contributions 
of higher-order multipoles. More significantly, the effects of 

external shear have not boon included for several of the sys- 
tems they analyze. We now set out to address these two 
concerns. 

In the case of arbitrary fcmax, SVD produces a particular 
solution x'"^ as well as a basis for the null space of A: {f }. 
We then have a family of solutions. 



A') 



(21) 



where Np = (3n — 1) > Nc = (2fcmax + 1) are the numbers of 
parameters and constraints, respectively. We must now se- 
lect appropriate coefficients d in order to construct the most 
plausible solution. Since the lens galaxies we are considering 
are elliptical, it seems reasonable to find the model with the 
smallest deviation from elliptical symmetry. In other words, 
wc want to minimize the wiggles in the isodensity contours. 

For a curve of constant k, the deviation 6r{9) from per- 
fect elliptical symmetry is given by: 



6r{9) 



where 



r{e)-ro{9), 

^ fcmax 

— \ (1 — (ttfc cos kO + bk sin kO) . 

fc=3 



(22) 



y - 3(02 cos 261 + 62 sin 261) 



is an isodensity curve for a perfectly elliptical galaxy. To 
quantify the wiggles, we average Sr^ over 9: 



Sr., 



{Sr')e = -^Y.^l-ky{al+bl). 



(23) 



fc=3 



Wc arc interested in the solution for which the root mean 
square wiggle is minimized. 

If we consider higher order multipoles but ignore shear, 
wc simply need to minimize the RMS wiggle with respect to 
the coefficients Ci. Since (Sr^) is quadratic in Ofe and bk, and 
hence also in a, this minimization is straightforward. 

When shear is included, the task of minimizing {Sr^) 
becomes slightly more involved. In particular, the lens po- 
tential of equation (3) must be modified: 



ip{r,9) = rF{9) - ^r"^ cos 29 - ^r^sin26>. 



(24) 



for constants 71, 72- This modification requires that the func- 
tions 7fc and Sk of equation (15), and the vector b of equation 
(18) be redefined by the expressions of Appendix D of Evans 
& Witt (2003). Namely, 



-TkiOi) 



(1-fc^) /„in cos fc6'nW^(6>„) 



Ski9i) = (1 - fc') 



-r„ cos fc6ii 1^(61;; 

fniri sin k9nW{9n) 
-r„ sin k9iW{9i) 



(25) 



(26) 
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where 

(6)0 = 1 + 71 cos 261; + 72 sin 201 (27) 
and 

b = [(1 + 7i)a;i +72yi, . . . , (1 + 7i)x„ +722/™, 
(1 - 7i)yi +72a;i, . . . , (1 - 7i)j/„ + 72X„, 
(/ni - l)rnri(l - 7? - 7I), . . . , 
(/„,„_! - l)r„r„_i(l - 7? - 72')]. (28) 

We see that including (71,72) as parameters to be deter- 
mined leads to a set of non-linear equations. To deal with 
this problem, wo use a non-linear optimization procedure. 
For specific values of (71,72), we can use SVD along with 
the minimum wiggle criterion to solve for the source posi- 
tion and Fourier coefficients. A minimization function can 
be employed to find the optimal values for (71, 72)- We refer 
to the resulting solution as the minimum wiggle model. 



3 RESULTS 

Let us now apply our methods to five quadruply imaged 
quasars. We begin with the Einstein cross, Q2237-I-0305, 
which docs not exhibit anomalous flux ratios at radio wave- 
lengths, thus providing a simple test case for the multipole 
expansion approach. Figure 1 shows a model with femax = 5 
and no external shear (cf. Evans & Witt 2003, Figure 2). 
The model exactly fits the observational data presented by 
Falco et al. (1996) and the CASTLES website. The eUipti- 
cal appearance of the isodensity contour confirms that the 
multipole method can find reeisonable solutions in lens sys- 
tems that do not require small-scale structure. Now let us 
turn to four anomalous systems: B1422-I-231, B2045-I-265, 
B0712-I-472, and B1933-h503. 

3.1 External Shear 

To motivate the need for external shear, we first study 
B1422-I-231 using a model that does not include shear. To fa- 
cilitate comparison with the results of Evans & Witt (2003) , 
we consider a fifth order multipole model, which fits the data 
of Impcy ct al. (1996) and Patnaik et al. (1999) exactly. It is 
clear from Figure 2(a) that this model is unphysical. In ad- 
dition to having a completely nonelliptical appearance, the 
model isodensity contour crosses the origin indicating that 
r becomes negative. Evans & Witt (2003) found a slightly 
better solution by inflating the measurement uncertainties 
(see their Figure 6). Even so, the model is incompatible with 
the observed structure of elliptical galaxies, a point noted by 
the authors. 

When we extend the model by including shear, it be- 
comes possible to find a physical — albeit not necessarily 
plausible — angular structure for the lens galaxy (see Fig- 
ure 2(b)). Before we can draw any conclusions from this 
result, however, we must determine whether the shear pa- 
rameters and Fourier coefficients we obtain are compatible 
with other observations. The shear parameters we find for 
B1422+231 are reasonable, because the lens lies in a group 
of galaxies that create a strong tidal field (see Kundic et 
al. 1997; Momcheva et al. 2005). In particular, our shear 
amplitude of 7 s •\/7i + 7I = 0.172, and orientation 



1 




• A 
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Figure 1. Isodensity contour with k = 1/2 (coincident with criti- 
cal curve) for Q2237-I-0305 with femax = 5 and no shear. The axes 
are labeled in arcseconds. 



0^ = (l/2)tan-i (72/71) = 40° are similar to those quoted 
by Kundic et al. (1997). 

The Fourier coefficients through order 4, plus some 
other model properties, are given in Table 1. To interpret 
them, we can determine the dimensionless octopole ampli- 
tude, A4, which describes the boxiness or diskiness of the 
isodensity contours, and compare it with the octopole ampli- 
tudes measured for the isophotes of elliptical galax;ies. The 
comparison is not perfect because lens models involve the 
mass while observations involve the light, but we can at 
least get a sense of whether the lens models are reasonable. 
The octopole amplitude A4 is just the Fourier coefficient for 
the density, expressed in a coordinate frame aligned with the 
major axis of the galaxy, and normalized by the semi-major 
axis length. In terms of the coefficients in Table 1, the major 
axis lies along the angle 62 = (1/2) tan^^(62/a2). Rotating 
into this coordinate frame then yields 

15 

A4 = - — — (a4 cos 46I2 + 64 sin 46I2) . (29) 

(The factor of —15 comes from 1 — , which appears when 
we convert from Fourier coefficients in the potential to those 
in the density.) If A4 is negative (positive), the isodensity 
contours are boxy (disky). For B1422-f231, our model with 
fcmax ~ 5 has A4 — —0.056. For comparison, the octopole 
amplitudes for the isophotes of real elliptical galaxies are in 
the range -0.015 < A4 < 0.045 (Bender et al. 1989; Saglia 
et al. 1993). In other words, the fcmax = 5 model is much 
more boxy than real galaxies (which is not surprising in 
light of Figure 2(b)). If we increase femax to 9 (see below), 
we find A4 = —0.021 which is still rather boxy. Going to 
femax = 25 yields A4 = —0.010, which is no more boxy than 
many observed elliptical galaxies. 
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(b) Nonzero shear 



Figure 2. Isodcnsity contours with k = 1/2 (soUd) and critical curves (daslicd) for B1422+231. Panel (a) shows a model with fcmax = 5 
and no shear. Panel (b) shows the solution for the same value of fcmax, but nonzero shear parameters (71,72) = (0.029,0.170). The dots 
show the image positions (suppressed in panel (a) for clarity). Note that with nonzero shear, isodensity contours and critical curves are 
not identical. 
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Figure 3. Isodensity contours (solid) and critical curves (dashed) for minimum wiggle models of B1422+231. Panel (a) shows a model 
for fcmax = 9 and shear parameters (71,72) = (0.028, 0.175). Panel (b) shows a model for femax = 25 and shear parameters (71,72) = 
(0.030, 0.167). 
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Table 1. Fourier coefficients, normalized RMS wiggle, and shear parameters for the various multipole lens models discussed in this 
paper. We quote models that fit the data exactly. 



System 1 




B1422-H231 




B2046-I-26B 




B0712-H472 




1 B1933-I-503 


fcmax 


6 


9 


25 


9 


17 


5 


6 


7 


8 


flEin 


0.797 


0.781 


0.779 


1.11 


1.11 


0.715 


0.710 


0.706 


0.515 


12 


-0.021:1 


-0.0159 


-0.0136 


0.00158 


-0.0211 


-0.0306 


-0.0424 


-0.0361 


-0.0112 


b2 


-0.0689 


-0.0492 


-0.0517 


0.0139 


-0.0135 


0.0394 


0.0246 


0.0239 


-0.0586 


as 


-0.00817 


-0.00359 


-0.00223 


-0.00928 


-0.000950 


-0.00333 


0.000379 


-0.0000154 


0.00144 


b3 


0.0132 


0.000983 


0.0000806 


-0.0159 


-0.00401 


-0.00595 


-0.000511 


0.0000860 


-0.00143 


04 


-0.00226 


-0.00119 


-0.000686 


0.000673 


0.000680 


-0.00466 


-0.000709 


-0.000234 


0.000869 


64 


0.00198 


0.000196 


-0.000196 


0.00367 


-0.000128 


-0.000546 


0.000113 


0.000287 


-0.000358 


<5rrms/ REin 


0.166 


0.0528 


0.0376 


0.235 


0.0943 


0.215 


0.0993 


0.0546 


0.126 


71 


0.029 


0.028 


0.030 


0.119 


0.092 


-0.044 


-0.072 


-0.066 


-0.033 


72 


0.170 


0.176 


0.167 


0.125 


0.090 


-0.078 


-0.092 


-0.085 


0.024 
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Figure 4. Isodensity contours (solid) and critical curves (dashed) for minimum wiggle models of B2045-t-265. Panel (a) shows the solution 
with femax = 9 and (71,72) = (0.119, 0.125). Panel (b) shows the solution with fe^ax = 17 and (71,72) = (0.092,0.090). 



3.2 Higher Order Multipoles 

Our next step is to include higher order multipole terms. 
In the case of B1422+231, the lowest order series for which 
a somewhat elUptical isodensity contour can be obtained is 
for fcmax = 9 (see Figure 3(a)). As we increase fcnmx the 
long-wavelength components of galactic structure disappear 
in favor of wiggles that are localized near the positions of 
the bright images A, B, and C (see Figure 3(b)). In other 
words, away from the images the model is smooth thanks to 
the minimum wiggle criterion. But there must be small-scale 
structure in the vicinity of the images in order to explain 
the observed flux ratios. Since the wiggles in the isodensity 
contour are not dramatic, it is not clear whether they should 
be interpreted as real features or just as approximations of 
other sorts of structure (such as CDM clumps). 

Let us now turn our attention to B2045-I-265. We find 



models that fit the data of Fassnacht et al. (1999) exactly. 
Unlike the case of B1422-I-231, the isodensity contour we 
obtain for fcmax = 9 is completely unreasonable (see Figure 
4(a)) and we must include multipoles of order 17 to ob- 
tain a remotely plausible model (see Figure 4(b)). Similar 
to B1422-I-231, we find that including higher order multi- 
poles reduces the RMS wiggle, but pronounced deviations 
from ellipticity remain primarily near the bright images. 
This suggests that the structure required to fit the anomaly 
in B2045-I-265 truly is local to the images. 

While the most obvious anomaly in B2045-I-265 is in the 
A/B/C triplet, it is worth noting that our models also have 
small-scale structure in the vicinity of the faint image D. 
Doblor & Keeton (2005) also concluded that the flux of im- 
age D is puzzling, and suggested that it has more to do with 
complex structure in the lens galaxy (such as an isophote 
twist) than with substructure per se. We cannot examine 



Multipole Models of Gravitational Lenses 7 




0.4 



1 1 1 1 
; • 


1 1 1 1 1 1 1 1 


1 i 1 1 1 
B1422H 


h231 


' 1 






















' 1 .... 1 .... 1 .... 1 .... 1 ' 


1 1 1 1 
_ • 




B2045H 


h265 


' 1 . 


• 


• • 








■ 1 . . . 










. 1 " 


5 


10 15 


20 




85 



B/A and C/A uncertainties to be 10% and 20% (we always 
use the observed flux uncertainty of 43% for image D). 

Since our formalism always produces models that fit the 
data exactly, the way we include measurement uncertainties 
is to add noise to the data and repeat our analysis. For every 
run, we perturb each data value by a random number drawn 
from a normal distribution with the appropriate dispersion, 
and then fit our model. We repeat this process 100 times, 
and select the case that has the smallest mean square wiggle. 
In this way wc find the miiiimurn wiggle model that fits the 
data within the measurement uncertainties. 

Figure 6 shows the results for B0712+472. As we would 
expect, including nicasurcmont uncertainties produces mod- 
els with slightly smaller wiggles. However, the changes are 
not significant enough to transform an implausible model 
into an acceptable solution. 



3.4 A Multi-Source Lens Model 



Figure 5. Normalized RMS wiggle {Srims / Reih) as a function of 



that hypothesis here because our models arc intrinsically 
self-similar, but it will be interesting to keep this image in 
mind as we develop more general models in the future. 

When examining isodcnsity contours, it may not bo 
completely obvious that high-order models with small-scale 
undulations really have a smaller total wiggle than low-order 
models. To understand that, recall that our minimum wiggle 
criterion (see eq. 21) is designed to select the model whose 
isodensity contours deviate least from elliptical symmetry. 
As femax increases, large-scale departures from ellipticity can 
be traded for smaller-scale features localized near the images 
in a way that does in fact decrease the total wiggle. Indeed, 
Figure 5 shows that Sr^ms/REin does decrease monotoni- 
cally as the order of the multipole expansion increases. The 
decrease is rapid at small fcmax but slows down as femax in- 
creases, and that gives us a sense of the order at which the 
multipoles have basically converged to have the minimum 
amount of small-scale structure. 



3.3 Measurement Uncertainties 

So far we have considered models that fit the data exactly, 
since that is the natural outcome of an SVD analysis of 
the underconstrained matrix equation (18). However, it is 
important to consider whether measurement uncertainties 
affect our conclusions. While we could do this analysis for 
B2045-I-265, we believe that the puzzling flux of image D 
would complicate the interpretation. We turn instead to 
B0712-F472. 

We assume the image positions to be precise to within 
±3 mas, which is slightly conservative compared to the 1 
mas uncertainties claimed by Jackson et al. (2000). For the 
flux ratios, we use the data of Koopmans et al. (2003), who 
found the uncertainties in the flux ratios of images B, C, and 
D relative to A to be 7.2%, 8.9%, and 43%, respectively. To 
be conservative, we construct models for which we take the 



To conclude this section, wc study B1933+503, where 
two sources arc lenscd by a single galaxy into two four- 
image conflgurations (Cohn et al. 2001, and references 
accompanying their Table 1). Both the image positions 
and flux ratios corresponding to one source arc known, 
while only the image positions of the second source have 
been determined to reasonable precision. We therefore have 
Nc = 2(ni -|- 712) -|- ni — 1 = 19 constraints, compared 
with Nc = 3n — 1 = 11 for the other systems we have 
analyzed, where rii = n2 = 4. Since wo must now flt a 
second source position, our list of parameters increases 
by two: Np = 2fcinax + 3. The matrix A then has the 
dimensions [2(ni + 712) + ni — 1] x (2fcmax + 3), and is given 
by adding more rows that represent the additional position 
constraints, as follows: 



A = 



10 

10 

10 



"01 



10 a, 



(2) 



10 

1 







'0,n, 



.,(1) 



"21 



(2) 
*2n2 



:,(2) 



(1) 



^(2) 
-^21 



o(2) 



5(1) 

2,n, 



(30) 



We first consider just the primary set of four images, 
which are labeled 1, 3, 4, and 6. (There is a flux ratio 
anomaly such that image 4 is brighter than expected.) The 
minimum wiggle model with fcmax = 8 that fits these images 
is not very plausible — it has a large protrusion near im- 
age 4, and a smaller one near image 1 — but at least the 
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(a.) femax — 5 



(b) femax = 6 




(c) /Smax — 7 



Figure 6. Model isodensity contours for B0712+472 with increasing multipole order. The sohd curve in each panel fits the data exactly. 
The short and long-dashed curves fit the data to within 10% and 20% flux uncertainties, respectively (except for faint image D, whose 
flux uncertainty is fixed at the observed value of 43%; see Koopmans et al. 2003). 



density is positive definite. When we add the additional con- 
straints from the positions of images 2a, 2b, 5, and 7, they 
dramatically reduce the solution space (fcmax = 8 is the low- 
est order case that is not over-constrained) . The only models 
that remain are unphysical. Modestly increasing fcmax docs 
not help. In other words, multipole models cannot simulta- 
neously fit the anomalous fluxes of images 1, 3, 4 and 6, and 
the positions of images 2a, 2b, 5 and 7. 



4 CONCLUSIONS 

We have shown that extending the method of Evans & Witt 
(2003) by including external shear and higher-order Fourier 
terms is essential for understanding whether multipole mod- 
els can fit observed lenses. The results for B1422-I-231 we 
have obtained are of a very different character from those 
of Evans & Witt (2003) . Even so, we cannot conclude that 
the multipole approach provides an acceptable explanation 
of flux ratio anomalies. 



The system B2045-I-265 requires multipoles of order 
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(a) Including image positions and flux ratios for images (b) Including all constraints 

1, 3, 4, and 6 only 



Figure 7. Isodensity contours (solid) and critical curves (dashed) for two models of B1933+503. Panel (a) includes the constraints for 
images 1, 3, 4 and 6. Panel (b) shows a model in which all of the observational constraints have been fit. Both models account for external 
shear and multipoles up to femax = 8. 



> 15. Even this level of small-scale structure leads to a 
rather wiggly angular dependence of density near the three 
bright images. Next, our analysis of B1933+503 reveals that 
a fundamental difficulty exists in fitting a multi-source lens 
with a simple multipole model. Finally, our method natu- 
rally finds models that fit the data exactly, but we have 
shown that our conclusions are not very sensitive to mear 
surement uncertainties. 

Our results suggest that there is a more fundamental 
problem with the global approach taken by Evans & Witt 
(2003) and ourselves in the current paper. It is possible 
that the problem simply comes from our choice of small 
scale structure. Sines and cosines provide a useful but by 
no means unique basis for carrying out a series expansion 
of the angular part of the potential. In addition, we have 
assumed a particular form for the radial dependence — no- 
tably self-similarity — which may need to be modified in 
order to find an acceptable galactic density function. While 
the present analysis seems to rule out simple multipole mod- 
els, the question of whether CDM clumps provide the only 
plausible solution has yet to be fully answered. 
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